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Abstract 

Wavelet bases are used to generate spaces of approximation for the resolution 
of bidimensional elliptic and parabolic problems. Under some specific hypothe- 
ses relating the properties of the wavelets to the order of the involved operators, 
it is shown that an approximate solution can be built. This approximation is 
then stable and converges towards the exact solution. It is designed such that 
fast algorithms involving biorthogonal multi resolution analyses can be used to 
resolve the corresponding numerical problems. 

Detailed algorithms are provided as well as the results of numerical tests on 
partial differential equations defined on the bidimensional torus. 
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Introduction 


Variational approximation methods for partial differential equations are based 
on weak formulations and on suitable spaces of approximation. 

Wavelets are known to be unconditional bases for a large variety of spaces 
and therefore are good candidates for the generation of approximation spaces 
for partial differential equation problems. The goal of this paper is to show that 
moreover, wavelet bases may lead to fast and adaptive numerical resolution of 
the corresponding approximations. 

In this paper, as in previous papers (J. Liandrat and P. Tchamitchian [10], 
[11]), the wavelets are used to expand the approximated solution of a partial 
differential equation as well as to approximate the kernel of the differential 
operator. They are not used only to perform “the linear algebra” (G. Beylkin 
[2]) related to more classical methods of resolution. 

Starting with an expansion of the form / = J2x < f^x > \£a, the solution 
of the equation Lu = / where L is a constant coefficient elliptic differential 
operator is u(x) = / f(y)K(x,y)dy where K(x,y) = J2x 

Under suitable conditions that will be made precise later, the functions 
as well as L^\ are pseudo wavelets, very close to wavelets (Y. Meyer 
[13]). This turns out to provide a stable approximation of u. However, the effi- 
ciency of the corresponding numerical approximation of u relies, at least in this 
work, on the hierarchic structure of multiresolution analysis since it provides 
fast tree algorithms. We will show that, if the operator satisfies suitable condi- 
tions that will be made explicit later, then the above mentioned pseudo wavelets 
are directly related to biorthogonal multiresolution and wavelets. Under these 
conditions, competitive numerical algorithms involving 0(N ) or 0(N log N) op- 
erations can then be derived. 

This paper provides the analysis of the problem and exhibits the correspond- 
ing numerical schemes. It is then organized as follows. 

The first part is devoted to the general concept of biorthogonal multiresolu- 
tion analysis on L 2 (IR n ). In the second part we focus on the problem of the sta- 
bility of the multiresolution framework under the action of constant coefficient 
elliptic operators. The cases of homogeneous and non homogeneous operators 
are treated separately. The third part deals with the numerical algorithms while 
the last section is devoted to numerical tests related to the resolution of elliptic 
and parabolic equations in bidimensional spaces. 


I Generalities: Biorthogonal Multi- Resolution 
Analysis in L 2 (JR n ) 

The concept of multiresolution is at the basis of our construction and we there- 
fore start with a short description of it : 
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Definition 1.1 (Y. Meyer, [12]) 

A r-regular multiresolution analysis of L 2 (IR n ) is a sequence of increasing 
closed subspaces Vj , j £ % , Vj C Vy+i, satisfying the following conditions 

i) ftz Vj = {0}, Vs is dense in L\JR n ); 

ii) f(x) £ Vj f(2x) £ Vj+i; 

Hi) f{x) £Vo<=> f(x - k) £ % Vfc € % n ; 

iv) there exists a function $ in Vo, such that the set of functions {$(#*-&), k £ 
Z n }, is a Riesz basis t for Vo; 

v) the function is regular and localized : $ is C r ~ l , is almost 

everywhere differentiable , and for almost every x £ IR n , for every integer 
a < r and for all integer p, it exists C p such that 

\d«<S>(x)\ <C p (l + J*j)^. (1) 

A consequence of ii), Hi) and iv) is that each Vj is generated by the family 
of functions {3>j*(x) = 2i n l 2 $(2lx — k), k £ Z 71 }. 

For simplicity reasons, we will only consider the case n = 2, but all the 
results presented to this article can be generalized in any dimension. We will 
always use for vectors a contracted notation: if e is a bidimensional vector then 
€ = (ci,e 2 ). 


1.1 Orthogonal multiresolution analysis 

To build an orthonormal multiresolution analysis, the Riesz basis {$(. — k), k £ 
Z 2 } is first orthonormalized in such a way that the resulting orthonormal basis 
is still of the form {$(. — k), k £ Z 2 }. 

The wavelets are introduced via the orthogonal complement of Vq in V \ : Wo . 
More precisely, if E is the set of all vertices in the unit cube [0, l] 2 , and if 
E* = E\ {0}, we have the following theorem: 


Theorem LI There are 3 functions W , e £ E* , in Wo, such that the collection 
{y$ e (x — k), k £ Z 2 ,e £ E*} is an orthonormal basis of Wo. Moreover , each 
satisfies the same property (1) of regularity and localization as $ and, moreover , 
satisfies the following cancelation property 


3 m £ IN , such that V k = (ki, fc 2 ) € 2Z 2 , 0 < k{ < m, 


[ x k V £ (x)dx = 0 . 

Jm? 


*A collection of vectors {e*, A € A}, in a Hilbert space H is a Riesz basis if any vector 
x £ H can be written in a unique way as x = ]P a A e A where (]P la *] 2 ) 1 / 2 is finite and 
defines a norm equivalent to ||n||jj. 
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The scaling invariance property ii) implies that, for all the family k E 
2£ 2 , e E E*} is an orthonormal basis of Wj. We will also use the following 
contracted notation: A E A^} where A j = {A = (& 4- 1) , k E ^ 2 ,£: E 

2?*}. Indeed, there is a straightforward bijection between A j and the set of pairs 
{(e, k ), k E E 2?*}. We will also use the following sets: A = Uj^zAj and 

A«=UjL n A i . 

From the inclusion Vo C Vi the following scaling relation can be derived: 

■*(*) = E h ‘*( 2x - o > 

J&Z 2 

while from Wo C V\ one obtains the following detail relation : 

&(x) = 9i *{ - 0. Ve € E* . 

IGZ 2 

It is very useful to transform these relations using the Fourier transform 
which is given by the equality 

/(O = f }{x)e~ i ^dx. 

Jm 

Indeed, the scaling relations then become 

®(0 = Mo(t/2)kt/2) (3) 

and 

r (^) = M e (£/2)%/2) (4) 

where Mq(£) = YLx^z 2 hie"^' 1 ) and M e (f) = YjI^z 2 flf are C 00 2 tt peri- 

odic functions. 

This leads to : 

_ 00 c _ c 00 c 

m = n^) , and no = (5) 

i=i i= l 

The following conditions are satisfied 

£ **«(£ + ™)MAZ + «) = ««' , V (c,0 € £ 2 (6) 

and 

Me(7re) = 4i,e 1 5 e J ,e 2) V(.£.,e) € £ 2 (7) 

and are called, following electrical engineering terms, the quadrature mirror filter 
conditions. We will also call the functions M e quadrature mirror filters. 
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Conversely, it is shown in A. Cohen et al. [4] that four 2i r periodic functions 
M e (£), 6 G & satisfying the quadrature mirror filter conditions (6) generate 
through (5) an orthogonal multiresolution analysis if some specific conditions 
are satisfied. 

Remark: 

• In this paper, we will often use specific multi-resolution analyses of L 2 (M 2 ) 
based on a tensorial product of multiresolution analyses of L 2 (1R ). More 
precisely, such analyses are defined as follows: if (Vj) is the sequence of 
spaces of a ID multiresolution analysis and if Wj, ip, ra 0 and m\ are re- 
spectively the related wavelet spaces, the scaling function, the associated 
wavelet and the quadrature mirror filters, then, the sequence of spaces 
( Vj ), defined as Vj ■== Vj 0 Vj is a multiresolution analysis in 1R 2 . More- 
over, <3>(a:) = <p(xi)<p(x 2 ) is the corresponding scaling function; (Wj), with 
Wj = W j are wave l e t spaces; the three generating wavelets are 

*o,i(z) = <p(x x)rp(x^), ¥i,i(s) - ip(xi)ip{x 2 ) and 9i,q(x) = ip(x 1 )<p(x 2 )] 
M £ (£) = rn £l (^i)m e2 (^ 2 ) with e € E are the quadrature mirror filters. 

L2 Biorthogonal Approach 

A relaxation of some properties of orthogonal multiresolution analysis can be 
performed using the one biorthogonal approach. This approach provides some 
flexibility since it allows to distribute the relevant properties of the multires- 
olution (number of zero moments, compact support or regularity) to the two 
involved multiresolution analyses. Moreover, it will turn out to be that the 
biorthogonal framework is “stable” under the action of a large class of opera- 
tors while the orthogonal framework is “fragile”. 

Definition L2 We call biorthogonal multiresolution analysis of L 2 (IR 2 ), two 
multiresolution analysis ( Uj)j£z and (Uj)jez such that there exists two families 
of corresponding scaling functions r and r such that: (Tj ki Tj‘ k l ) = SjjtSkk' for 
all j, f £ % and , k and k ' £ E 2 . 

In this case we define the wavelets spaces Xj and Xj as Uj+\ — Uj ® Xj , 
Uj+i = Uj ® Xj with UjLXj, UjLXj , and we introduce the functions 9\ — 
2^9* (23 . — k) and 6\ = 2 J 0 £ (2 J . — &), s £ E , that generate respectively Xj and 
Xj and such that {#Aj 9\t) = 

Moreover, following the construction of orthogonal multiresolution analysis 
we define 2x4 filters (i.e, C°° 2tt periodic functions), P £ and P £ , e £ E ) 
associated with the two biorthogonal multiresolution analyses. These filters 
satisfy the biorthogonal quadrature mirror filter relations equivalent to (6) that 
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are 


\/e ,e' and e' € E, £ € [0, 2tt] 2 , 

p *(f + 7 re)P e /(£ + ?re) = <5 ee < , (g) 

€€£; , _ 

P £ (^Tte / ) — ^e.i',ci ^£2,621 ■Rr(^C / ) 3=1 ®€-i-,fii^ea,ca * 

As in the orthonormal case, the generalizations of relations (3) and (4) relate 
the scaling functions and the wavelets to the filters as 

oo 

m = Poiamsm = n p °( 2-i o, no = m/mm, (9) 

j= i 

_ ^ OO ^ 

m = P 0 (£/2)?(£/2) = n p o(2" J '£), 9 (£) = P E (£/2)r(£/2). (10) 

j = l 

As in section I.l,Jhe question to know is under which conditions the biorthog- 
onal filters P £ and P € satisfying the quadrature mirror filter conditions (8) de- 
fine two biorthogonal multiresolution analyses? Again, a specific condition is 
required and has been formulated in A. Cohen et al. [4]. We will use a weaker 
version of this formulation adapted to the case of functions with fast decay. It 
can be expressed in term of the following theorem for which a complete proof 
can be found in Pj. Ponenti [16] and W. Dahmen and A. Michelli [5]. 

Theorem 1.2 Let a > 0 and let P e (f) and P £ (t ;), e £ E } be eight C 01 2tt periodic 
functions satisfying the biorthogonal quadrature mirror filters conditions (8). 
Defining r, r, 0 , and 6 using formula (9) and (10), if 

- there exist C and fi > 0 such that for all £ £ IB? 


m\ < c( i+m- 1 -*, 

m\ 

< c(i+ leir 1 -" , 

(11) 

- V k ,et k' , 




{ t(x — k),r(x 

-k')) 

= hk', 

(12) 

- and if 




JjR> l^l 2 ^ ’ 

s 

Jm 2 

fe)l 2 

|£|2 ^<°° > 

(13) 


then 


• the sequences of subspaces (Uj)jzz, and (Uj)j^z generated respectively 
by {Tjk,k £ 2h 2 } and {rj^^k £ 2Z 2 } are two biorthogonal multiresolution 
analyses; 

• the wavelet families { 0\{x) == 2^ £ (2^x - k), A £ A } and {0\(x) = 
2?6 s {2?x — Ar), A £ A} are two biorthogonal Riesz bases of L?{1B?). 
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II Constant Coefficient Elliptic Differential Op- 
erators and Biorthogonal Multiresolution Anal- 
ysis 

II. 1 General results 

The starting point of this section is the following remark. Given (\Pa) a family 
of orthonormal wavelets and knowing / = a)^a> the solution of the 

equation 

Lu = f, (14) 

where X is an elliptic operator of order s is, at least formally, 

u = J2(f^4 L ~ l l^] ■ (15) 

When X and X” 1 are bounded on X 2 (1R 2 ), the families {L~ ll $ x } and {X*W a} 
are two biorthogonal Riesz bases. 

The question we address now is related to what happens in the specific case 
of wavelets when the operator X is unbounded (as in the case of a differential 
operator). 

In the following paragraphs, we first see that, assuming some compatibil- 
ity conditions between and X, the two families of functions {L“ 1 $a} and 
{ZA^a} are wavelets or pseudo wavelets (Y. Meyer [13]) depending on whether 
the operator is homogeneous or not. Then, we show that in some cases, a 
biorthogonal framework embedding {X” 1 ^} and {L *\ If a } can be built. 

To be more precise, we take (Vj) an r-regular multiresolution of L 2 (M 2 ) con- 
structed using a tensorial product of two ID multiresolutions and X a constant 
coefficient elliptic differential operator. Let us write X = J2o<a<s a aD C( where 
D is the operator We define in a standard way the symbol of X as 

( 16 ) 

It is a polynomial in £ of degree s and, if / e C°°(1R 2 ), we have the well known 
formula 

V£ € JR 2 ,i/(() = /(M)- (17) 

We note formally 

9 £ = X*$ £ and 0 £ = X" 1 ^ (18) 

and more generally, 

0 £ x = 2“ J5 X*^a and 0\ = 2 js L^ e x . (19) 

Note that when X is homogeneous, 9 x (x) = 2^0 £ (2^x — k) and 9 x (x) — 

23 0 s (2,3 x — k) while, in general this is not true. 

Then we have 
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Theorem IL1 Given a family of r regular wavelets with m + 1 zero moments , 
if L is a homogeneous operator or if L is an inhomogeneous operator with a 
strictly positive symbol of order j > 0, and such that r > s and m > s, then 
Regularity: 0 e E ff r - s and Q £ £ H r+S 

Localization : for all multi-indices* and 7 such that \j r \ < r — s and 
M < r + 5 ~ 1 and integers l E IN , 

( |a 7 '^(x)| < CV 2 Jh ' l 2 J (l + 2-'>-A|) - ' 

\ |a 7 ^(x)| < C y 2^21 (l + 2f|x-A|)" 3+5 " ro_l71 

Cancelation: Let x k — x \ 2 ), 

I x k O e (x)dx = 0, 0 < |fc| < m + a 

i 

I x k 0 s (x)dx = 0, 0 < \k\ < m — a , 

Jx n 2 

where a — $ in the homogeneous case and a == 0 in the inhomogeneous one . 

A complete proof of this theorem can be found in Y. Meyer [13] and Ph. 
Tchamit chian [17]. 

Remark: 

• Following Y. Meyer [13] and according to Theorem ILl, a factorization 
of the operators L* and L~ 1 can be performed as L* = C o T~ s and 
L~ x ■= C g P where T is a diagonal operator in the wavelet basis defined 
as: ipx t — 2 J ij>\ and where C and C are bounded on L 2 and defined by 

C : 6\ , C : Ox , (20) 

The operators C and C act just as a transformation between two bases. 
The operator F is nothing else but the classical preconditioning operator 
for elliptic problems (S. Jaffard [8]) that mimics a diagonal derivation in 
the wavelet basis. 

In other words, thanks to this factorization, the computation of the image 
of a function by an elliptic operator or its inverse can be transformed into 
a well conditioned problem using a diagonal operator in a suitable wavelet 
basis. This is essential since it provides the numerical stability of the 
further developed algorithm. 


We can then rewrite (15) as 

u = Y l 2~ js (f, *x)9 x . ( 21 ) 

t We call a multi-index any couple of integer 7' = (7i,72)and d " 1 ' 6 = 
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An important issue, as far as numerical applications are concerned, is the 
computation of the sum (21). Indeed, even if this sum corresponds to a pseudo 
wavelet decomposition, fast algorithms for the computation of the sum are not 
available. In the framework of this paper, the fast algorithms are linked to 
the concept of multiresolution analysis presented in section I. We prove in 
the following sections that, under suitable conditions, the construction of a 
multiresolution analysis embedding the function 0\ is possible. Moreover, we 
provide explicit expressions of the quadrature mirror filters required for the fast 
implementation of (21). 

The starting point of our construction is due to P.G. Lemarie [9] who con- 
structed in the one dimensional case two biorthogonal multiresolution analyses 
from an original orthogonal one and from the derivata operator. We generalize 
this approach to any dimensions and for any homogeneous elliptic operator. 

In contrast to the classical constructions of multiresolution analysis, this is 
an inverse problem. Indeed, knowing the two dual wavelet bases we can define 
two sequences of subspaces (Xj) and (Xj): 

Xj = span{f?A , A € Aj , ./ < j} ^ 

Xj = span{0A, A gAi,/ < j} , ' ’ 

The open question is the following: how can we construct two sequences of sub- 
spaces (Uj) and (Uj) for which (Xj) and (Xj) play the role of two biorthogonal 
wavelet spaces? 

In other words, the problem is the construction of the generalized scaling 
functions r and r related to 0 and 0. 

In the case of non homogeneous elliptic operators the approach used in the 
homogeneous case can not be transposed and we will not be able to define a 
multiresolution framework embedding the space sequences (Uj) and (Uj). How- 
ever, we will show that the essential property of embedding spaces as well as 
the existence of scaling (3) and detail (4) relations can be saved. This will allow 
us to derive fast and stable algorithms to sum up (21) even in the case of non 
homogeneous operators. 

11*2 The case of homogeneous elliptic operator 

In that case the natural candidates for r: Zr 1 ^, are not defined in L 2 for the 
basic reason that $ does not belong to the range of L, or, in other words, suffers 
from a lack of zero moments. Using a preconditioning operator, we will adapt 
the function <3> to the operator L~ l (i.e. we will transform $ so that it enters 
the range of L) while preserving the two scale relations (3) and (4). We will 
then check that the resulting multiresolution analysis is fitted to the functions 
9\ and 0\ previously defined. 

More precisely, we have the following theorem: 
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Theorem II. 2 Let (Vj) be the family of embedded spaces of an r -regular mul- 
tiresolution analysis ofL 2 (IR 2 ), let L be a homogeneous operator of order s and 
of symbol a and let S be a periodic function, not vanishing on ]0,27t[ 2 and 
equivalent to a in zero . Ifm>s + 1, r > 5+1, and if the eight functions P s (£) 
and P S (Q> e E E defined as 

MO = 2 s == Mo(0, Po(0 = L Mo(0 , (23) 

mo = 2 Mm(o, mo = 57^ (24) 

are C a , a > 0 then they satisfy the quadrature mirror filter conditions (8) and 
define two biorthogonaljnultiresolution analyses. The corresponding wavelets 
are the functions 9 and 9 and the scaling functions r and r are derived following 
(9) and (10). 

Proof: 

□ By construction all the filters are C a with some a > 0 and they satisfy 
the biorthogonal quadrature mirror filter conditions (8). 

The only point to prove is the convergence in L 2 {1R 2 ) of the infinite prod- 
ucts (9) and (10) defining the two scaling functions. We will use the following 
lemma: 

Lemma II.l Let p(x), x E M 2 , be a homogeneous polynomial of degree s, 
and let S and C, be 2tt periodic functions; then the following propositions are 
equivalent: 



J5 cr '' ) = X) ■ 

(25) 

a) 

S(x) = 2 s C(x/2)S{x/2) , 

(26) 


S(x) ~ p(x) . 

X — ►u 

(27) 

Proof: 

□ The equation (26) is obtained from (25) written for x and x/2. 7 while (27) 
is derived from (25) when x — ► 0 since necessarily C(0) .== 1. 

Conversely, (25) is obtained from (26) . Indeed, since 


S(x) S{x/2) S(2- N x) jr • 

p ( ») - C( * /2 W) - Pi 2-"«) jy[ c(J *>• 
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thanks to (27) we obtain (25) when N -+ oo. ■ 


This lemma allows us to calculate the infinite product (9) and (10), and to 
get: 

oo oo oo 

m = n 

j = 1 

OO OO OO 

m = n p °( 2 ~ j t) = n^n^)«) , 

i= i i=i i=i 

and finally, 

? ( o = ||| m , H m • ( 28 ) 

The function $ being r-regular the conditions (11) are trivially satisfied. 
Defining the function 0 et 0 by (9) and (10) we get 

H) = P*W 2 )^/ 2 ) = ^( 0 %) 

~ s ~ $(£) 

e (o = p&mm = H • 

The wavelet admissible condition (13) is immediately satisfied thanks to Theo- 
rem II. 1. Finally the assumption (12) is satisfied by construction, that completes 
the proof. ■ 

Remarks: 

• Note that here, thanks to homogeneity, the subscript a recovers its classical 

“wavelet” meaning since in that case, 6\(x) = x — k) and 0\(x) ■= 

yl 2 e{{Tx - k). 

• The relation (25) is a generalization in any dimension, of the classical 
formula 

freest *) = 

r . x 

j-i 

used by P.G. Lemarie for the first order differential operator in one di- 
mension. 

• The function S can be interpreted as the symbol of a difference operator 

that we will call D l- If S is a trigonometric polynomial, then Dl is a 
finite-difference operator and S is C°° . The fact that S(x ) <j{x) is 

just the translation that Dl is consistent with L. Moreover, S removes 
exactly the singularity of r for £ == 0. Conversely for r, the singularity 
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given by S at points (2?m, 27m), n E will be removed exactly by cr in 
zero and by the zeros of $ at points (27m, 27m), n E Notice that this 
last point won’t be true for inhomogeneous operators. 

From a certain point of view, Dl can be seen as a preconditioner for L 
since r and r are defined by 

r = L~ 1 Dl$ and f = (30) 

• In one dimension there is a canonical choice for 5 and therefore for Dl 
such that, if the function <£ and have a compact support, then r, r, 
9 , and 9 are also compactly supported. Indeed, in that case we have 
necessarily a(£) :== a£ s , a E <P. Therefore, the canonical choice for S 
is £(£) — a(—i) s — and Dl is then a non-centered finite- 

difference approximation of L of order 1. Indeed, it is well known (see 
I. Daubechies [6]) that the quadrature mirror filters related to orthogonal 
compactly supported functions $ and $ are: mo(£) = (l -f- e^) £(£) 

and mi(£) = (l — e“^) m £(£ + tt) where C is a finite trigonometric 
polynomial. Then we get 

m) = 2 5 (1 + .e*)”’"' £(*), 

PliO = 2 s <mV' ? (1 - e-^) m (1 - e^) s X(e + tt), 

P 0 (O = 2~* (1 + e^) m (1 + £(£), 

?i(0 = 2-'.fl" (1 - e- £ ^) m S £(C + ») , 

which proves that P 0) -Pi, Po and Pi are also finite trigonometric polyno- 
mials. Then, using the following lemma borrowed from G. Deslauriers and 
S. Dubuc [7] we deduce that the functions r, r, 0, and 9 have compact 
support. 

Lemma II.2 

If m = £ with ZnlN^n = l then U.T=l T ( 2 ’ j O is «» 

entire function of exponential type. In particular, it is the Fourier trans- 
form of a distribution with support in [N\, iV 2 ]. 

Clearly, this canonical form is no longer available if the space dimension 
is larger than 1 since the multidimensional quadrature mirror filters can 
not be factorized as above. 

11,3 The case of inhomogeneous elliptic operator 

Here the non homogeneous property of the operator is obviously not adapted 
to the scale invariance property of the multiresolution analysis. We will see 
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however, that introducing at each level a new scaling function, an embedded 
family of spaces can be constructed which preserves the mathematical properties 
relevant for numerical applications. 

The natural candidates for r\ are a the functions They are well 

defined in L 2 (IR 2 ) but suffer now from a lack of localization when j increases. 
Indeed, we have 

. lira || T^L- X [<p jk ] (.*) - G(x - 2~ j k) || = 0 

j — -f+GO 
k 2”* J *Xq 

where G is the Greens function of the operator L defined as G(£) == l/er(£). For 
example, when L = 1 — A, G(£) = 1/(1 +£ 2 ) and G(x) = e"*^. G decreases fast 
but mathematical and numerical difficulties come from the fact that the family 
of functions {G(x — &2“ J ), k £ Z 2 } is not a set of functions rescaled with j (in 
other words, this family is not obtained by rescaling and translation of a single 
initial function). This implies that the control of the localization by the index 
j is lost. It follows that the family {L” 1 # a, k £ Z 2 } is not a good basis to 
reconstruct our solution. 

Let us show now that, however, a process very close to the one used in 
the homogeneous case will provide an efficient algorithm for the summation of 
formula (21). 

We mimic the construction performed in the homogeneous case. Let (Vj) 
be an r-regular multiresolution analysis, let L be an elliptic operator of order s 
with constant coefficients and a its symbol (we now suppose cr(£) > cr 0 > 0 V£). 
Let us also define the homogeneous polynomial of order $, <r, as the principal 
part of cr, and let 5(£) be a C°° 2tt periodic function with 5(£) Q € n , where 

n will be fixed later. 

Then, Vj £ Z, we define a difference operator Dj by its symbol S(£/2 ; ). 

Following the previous section, we define Vj , k £ Z 

r jk = 2 (31) 

and 


0 jk = 2 t'L-'&Jhh ( 32 ) 

By construction, and thanks to the fact that L is a constant coefficient 
operator we have Tjk(x) = Tj(x — ^) and Oj k (x) = 6j ( x — -—) where 


and 


= 2 , 
e s -(n = 

*(0 


(33) 

(34) 
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Remark 

♦ The functions fj mimic the function f defined in (28). Unfortunately, it is 
not possible to define the equivalent of r (28) since Df 1 L$jk £ L 2 . Note 
however that, by chance, (15) involves directly the 0\ functions. 

Then, with Po and P £ defined in (23) and 24 we get the following scaling 
and detail relations: 

|(0 = 2Po(t/2 j+1 jTj +1 (0 

0j(O = 2P £ (t/y+ 1 )T J+1 (t) , Ve € E*. 

Remark: 

• An important point is that the filters P £ are independent of the scale 

index j as it is originally the case for standard multiresolution analysis. 
Furthermore , since they are defined by (23) and 24, the filters P £ are 
directly related to the homogeneous operator L of symbol & if n — s. 
This point is essential since it means that if Dq is consistent to L, then 
the tree algorithms related to the multiresolution spaces (Uj ) and used to 
sum up (21) are stable even if the functions Tj{x) are not standard scaling 
functions. Indeed, S(£)/V(£) 0 and therefore r j (0) = 0. 

In other words, even if the functions Tjk(x) are used as scaling functions on 
the range of L, they have zero moments as wavelets have. 

Finally, we can prove the following theorem: 


Theorem II. 3 For 0 < n < ■$, s < r and s < m, the functions Tjk defined 
by (33) and (31) satisfy 

|fi rr r i *(*)l < c; (1 + 2 j \x - 2~ j k\y n+s ~ m ~ 3 ~ h[ . (36) 

If n = s and ifP £ (£) ; e E E are C a , then stable tree algorithms are available 
in QJ )£ := q -K-tj 

30 < C <C f < +oo such that if f = °^0a? then 

c£| 4 I 2 <||/|| 2 < c' £|dx|> . (37) 

A A 


Proof: 

• Since S £ C°° the we can apply Theorem II. 1 that proves the localization 
inequality (36). 
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• The relations (35) defines directly the tree algorithms required to compute 
the coefficients {cjk} such that 

^2 d x 9\ = 2J CjkTjk 
A6AJT 1 k 

and since, when n = s, the involved quadrature mirror filters could be 
related to the homogeneous operator L, the stability of the algorithm is 
equivalent to the stability of the transform {cjk} »—*■ / for / £ 

• This transform is stable if and only if the family {r^, k £ ?Z 2 } is a Riesz 
basis of 0^0 X t , 

We have 



Since {$jk} is a Riesz basis, since 5(£) is bounded, and since a is bounded 
from below, then || / || 2 < C r J2 \ c jk | 2 , which is the second part of the 
inequality (37). 

To prove the first part of (37) we use again the fact that the filters Po and 
P £ are related to L. 

Indeed, if we define 9\ replacing a by & in (29) and if we define / as 
f “ 2a then the transform / »-* {d\} is stable. Moreover thanks to 

theorem II. 1, the operator / >-+ f which can be also defined as *-*■ 6\ 
is also bounded (Y. Meyer [13]). Therefore the operator / {d\} is 

bounded that is the first part of (37) and that completes the proof. 


Ill Approximation and Numerical Resolution 
of Elliptic Problem on the Torus 

This section is devoted to the approximation of elliptic problems on a sequence 
of embedded Galerkin spaces associated with a multiresolution analysis and to 
the corresponding numerical algorithms. 
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Classically, we will consider the problem with periodic boundary conditions 
to avoid the difficulties of general boundary conditions. We will use a r-regular 
multiresolution analysis of the torus T[ o t ip = (M/2%) 2 constructed using a clas- 
sical periodization technique (Y. Meyer [12]). We take as granted that, with 
minor modifications, the results proved on the whole line can be transposed to 
the torus. In that section, homogeneous and inhomogeneous operators will be 
treated similarly. 

IIL1 General formulation 

The general formulation of the problem is 

Find u G 2|o,i] 2 such that 

Lu = f (38) 

with f € %ip and L a constant coefficient elliptic operator of 
order s. 

Standard variational approximation (P.A. Raviart and J.M. Thomas [14]) 
leads us to look for the solution of a weak problem in so called Galerkin approx- 
imation spaces V e , where e is a scale related to V e with ^ T[ 0) ip. 

A natural choice for V e is V e = Vf where Vf belongs to a multiresolution 
analysis of Tjo^p of the type described above. Indeed, we then have the following 
inequality guaranteed if the involved multiresolution analysis is r-regular, 

V s < r, 3 c > 0 V/ € H s | ]/ - U v vf\\ 2 < C \\f\\ n , (39) 

where Tlyr, j < 0 , stands for the orthogonal projection on V? . 

Then a standard Galerkin approximation writes 

Find u p G V* C T[ 0> ip such that 

^ ILyvLUyvUp = II yv f , (40) 

where II V v ) j < 0, stands for the extension operator from V? to T[ 0) ip . 

This approach leads us to replace L by the approximation UyvLUyv. The 
corresponding numerical algorithms are reduced to linear system solvers once a 
basis of Vp has been chosen (P.A. Raviart et al. [14]). 

§The symbol . v stands for the periodization operator on [0, l] 2 . We recall that dim Vf = 
2 2j = 1/3 dim Wj and that Vf* = span{# 0 o — 1}> 
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Numerically, an optimal choice for the expansion basis is the wavelet basis 
ofVJ: since the corresponding stiffness matrix, is sparse 

and can be easily (i.e using diagonal matrices) uniformally preconditioned (S. 
Jaffard [8]). 

111. 2 A different approximation 

The main purpose of this paper is to define a different approximation of L and 
the corresponding numerical algorithm. Taking the set of functions { 9 J , A G 
A^" 1 } (defined in (31)) and defining U p = span 0J, A £ A P }, we get 
an approximation of u as 

% = E + (/, • (41) 

A€Ag 

Indeed, 

u p = X- 1 nv P / = P^I- 1 / (42) 

where P~ v is the projection on U p orthogonal to U p . 

This formulation of u p shows that the* convergence of u p towards u when 
p oo is straightforward since the set U p is a family of Galerkin spaces for the 
suitable space of definition of u. 

Moreover, the stability of the algorithm is a direct consequence from the 
classical preconditioning properties of wavelet base expansions (S. Jaffard [8]). 

111. 3 General scheme 

The numerical algorithm derived from the previous section is now presented 
in its collocation version. We call Ij the set of points {7Z Q[0, 2 J [) 2 . Then, 
Jj = {2 ~ j k,k G Ij} is the two dimensional regular grid of scale 2~ J related to 

M 2 . 

For the numerical implementation, we assume that the space V p is such that 
any continuous function f G V p is unambiguously defined by its values on the set 
of points J p . This assumption (satisfied by the even order spline multiresolution 
we will use in the numerical tests) allows us to define the collocation projection, 
Cyv , from the set of real sequences {f k )kei p to V p as 

Wk)k€i p ,C v v((f k ) keIp ) = / &f€ Vp V and /( 2"'*) = /*, V* G 

As soon as we define / G V p by its coordinates on the basis ® pk , the opera- 
tor Cyv appears as a discrete convolution operator involving a so-called inter- 

p 

polant filter !$„(£)• The operator Cyl is also a convolution operator called the 

^ Again, AJJ, = U™_ n Aj with, in the periodized framework, Aj — {A = 2~ J (k+ |) ,k 6 

2Sf\ [°, 2 T) 2 ,ee£*}. 
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point value operator and involves the point value filter PV$ p (£) defined from 
{®po(Jp)}- Obviously, /$ p (£) and PV$ p (£) are inverse. Let us remark however, 
that the point value operator always exists as soon as the functions are 
continuous. 


For the implementation, we therefore replace IIy7> by C v v in (42) and define 
therefore u v as 

Up = lT l Cyvf . 

Given the point values of / on the grid points J p , the algorithm provides the 
values of u p on the same grid. More precisely, the algorithm can be presented 
as follows: 

1. The input of the procedure is the set of values (f(J p )) from which the 
interpolant function f p E V p is constructed using I$ p (£): 

f P = £ tp&h 

kei P 

2. f p is then decomposed into the wavelet subspaces Wf , 0 < j < p — 1 and 
V? as 

/ = £ (f,nm+coo 

A€A p 1 

where c 00 = U v v(f). 

3. u p then becomes 

Up = ^ + c oo 

AeA^ 1 

where 9% — (2i s L~ 1/ 4>\) v . Here, c' 00 ■= coo/cr(0) for non-homogeneous 
operators. For homogeneous operator Cg 0 should be given. Note that in 
that case cr(0) = 0 and / should have at least $ vanishing moments; the 
fact that Cqq should be given corresponds to the ill posed property of the 
initial problem in L 2 . 

4. u p is then expanded in terms of the set of functions {f^ ki k E I p } using 
the tree algorithms related to and Tpk as 

U P = 53 C P k ^Pk ‘ 

kelp 

5. Finally, the grid point values of u p on J p are estimated using the point 
value filter PV~ (£). 

T p 
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It appears clearly that various precalculations should be performed. In the 
first step, the interpolant filter related to must be known; for the second 
step, the orthogonal multiresolution analysis quadrature mirror filters M e should 
be used and, for the fourthstep the corresponding biorthogonal multiresolution 
quadrature mirror filters P s are required; finally, the point value filter related 
to r p is used for the last step. 

To be more precise, we have to make some remarks that help to reduce the 
complexity and storage. For steps one and two, tensorial properties can be used 
in a very classical way to reduce the 2D- algorithms involved in 2(dim(t^ p )) 1 / 2 x 
ID- algorithms. It is then enough to know the one- dimensional interpolant filter 
related to <fp k and the one dimensional quadrature mirror filters £ 6 E . 

For steps four and five, where the filters P €: and the point value filter PV~ are 
involved, we can note apply this simplification and we have a full 2D-problem. 

We are now able to summarize all these precalculations in the following step 

0 : 

0. The computation of the following filters is performed (this is presented for 
the spline multiresolution analysis case): 

-Interpolant filter related to ip p k , analytical formulas in one di- 
mension are available in V. Perrier and C. Basdevant [15]. 

-Orthogonal ID multiresolution analysis filters m € : analytical formu- 
las are also available. 

-Filters P £ of the Biorthogonal Multi resolution Analysis: These filters 
are constructed from m e and formula (23) and (24). In fact only jP(o,i) 
and P( i }1 ) have to be computed since we have P(i,o)(£i 5 £2) = ^(o,i)(£2>€i)- 

-Point value filter, PV — related to i This filter is computed from 

Tpk * 

formula (28) and the analytical expression of 5, cr and <£(£). We have 
successively 


?;(*) = 1/(2*) £ %(w)e 2i ™ * , x E R 2 , 
w£Z 2 

r;{z n ) = 1/(2*) £ ( Y, + 2Pr ) ) e 2i ™- Xn > *« € J P . (43) 

wEl p \r£Z 2 ) 

Practically, 43 is truncated according to a prescribed precision. This is 
possible because r p (f) decreases fast. 

Remarks: 

• One should again emphasize that the entire algorithm is based on con- 
volution operators. Thanks to the periodic boundary conditions, the 
convolutions can either be performed directly or using a discrete Fourier 
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transform. The implementation presented in this paper uses the Fourier 
transforms since it is optimal for non compactly supported filters on non 
adapted spaces of approximation. 

III. 4 Detailed Algorithm 

This section is devoted to the structure of the code. Basic tools, such as Fast 
Fourier Transforms, Convolution/ Decimation algorithms, or Term by term mul- 
tiplications are not described. 

As can be seen from the general scheme presented above, the main code 
involves only two more elaborate routines that will be called the Precalculus 
routine (step 0), and the tree algorithm routines. The tree algorithm routines 
mayor may not use the tensorial structure. They will be called consequently 2D 
Tensorial Tree Algorithm-D (steps 2) and 2D Non Tensorial Tree Algorithm-I 
(step 4) where -D and -I stand for direct and inverse. We recall that the steps 
1 and 5 are convolutions and the renormalization performed in step 3 is term 
by term multiplication. 

The tree algorithm routines are becoming very classical and therefore we 
will not describe them either. Note however, that since only convolutions are 
performed in our algorithm, we only use the discrete Fourier transform of the 
wavelet coefficients (and not the corresponding values) at every scale, that re- 
duces significantly the complexity. 

We now give the detailed description of the main program (ELLIP ) and of 
the precalculus program (PRECAL) in pseudo code. 

The following example sketches the structure our programs. 

[OUTPUTS] =Pr ograxa ( INPUTS ) 

# Comments 

# Body of Program: 

[OUTPUTS 1]- Subprogram! (INPUTS) 

INPUTS2 = OUTPUTS 1 
[OUTPUTS] = Subpr ogram2 ( INPUTS2 ) 

> 


Our variable descriptors bears some resemblance to the C language as well 
as to the MAT LAB conventions. 


III.4.1 Preliminary computations 

The symbols *, and ./ used to present this program are borrowed from 
MATLAB and mean respectively, the transposition, the matrix product, the 
term by term product, and the term by term division. We also use the following 
n sub sampling operator a : n : b defined as: If a is a 2D array of size 2 P x 2 P , 
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b(l : 2 n : 2?, 1 : 2 n : 2^) is a new array of size 2 p ~ n x 2 p "" n given by a(i 7 j) — 
a(2 n i, 2Nj), (i,j) £ l p . 

Program PRECAL 

[QMFBIW , TAUTW] = PRECAL (p , pmax , QMFW , PHIW , SW , CW , SIGMA ) 

#HfPUI : 

#p -> index of the approximation space Vp in which the 

# elliptic problem is solve. 

#pmax “> index of the approximation space Vpmax in which the 

# precomputation of TAUTILDE is done (it depends on the 

# prescribed precision) . 

#QMFW -> structure containing the quadrature mirror filters in 

# one dimension: 

# QMFW.mO “> ID array containing the quadrature mirror filter 

# coefficients associated to the scaling functions; 

# size ( QMFW . mO ) ->2"p ; QMFW.mO(i) = m0(i/2"p), 

# i belong to {0 , . . . ,2"p-l> . 

# QMFW.ml “> ID array containing the quadrature mirror filter 

# coefficients associeted to the wavelet; 

# size (QMFU .ml ) ->2"p ; 

# QMFU.ml(i) = ml(i/2"p) ,i belong to {0, . . . ,2"p-l}. 
#PHI¥ -> ID array; size (PHIW)->2 "pmax; where pmax is given 

# and pmax>p; PHIW(i) = the value of the Fourier transform 

# of the ID scaling function at the point i, i belong to 

# {0, . . . ,2"pmax-l}. Used to compute the value of tautilde on 

# the finer grid. 

#SU -> 2D array containing the sampling of the function 

# S used for biorthogonal filters; 

# Size(SW)->(2"pmax X 2"pmax) ; SW(i, j )=S(i/2"pmax, j/2"pmax) , 

# (i,j) belong to {0 , . . . ,2"pmax-l}"2. 

# 

#CW -> 2D array containing the sampling of the function 

# S(2w)/(2"s S(w) ) Size(CW) -> (2"p X 2"p) ; CU(i,j) = 

# S(2i/2"p,2j/2"p)/ S(i/2"p, j/2"p) , 

# (i, j) belong to {0, . . . ,2"p-l}"2. 

#SIGMA -> 2D array containing the sampling of the symbol 

# of the operator Size(SIGMA) -> (2 "pmax X 2 "pmax) ; 

# SIGMA(i,j) = sigma(i/2"pmax, j/2"p) , 

# (i, j) belong to {0, . . . ,2"pmax-l}"2. 

# 

fOUTPUT: 

#QMFBIW -> structure containing the biorthogonal quadrature 

# mirror filters related to the tautilde functions: 

# QMFBIW . PTILDE0-> 2D tab containing the biorthogonal 
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# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# 

# Computation of the filters Ptilde 

# 

QMFBIW.PTILDEO = ((QMFM.mO) > * (QMFW.mO)) .* CW; 

QMFBIW. PTILDE 1 = ((QMFW.mO) ’ * (QMFW.ml)) ./ 

(2~s SW(1:2~ (p-pmax) : 2~pmax, 1 :2~ (p-pmax) :2~pmax); 
QMFBIW. PTILDE3 = ((QMFW.ml) > * (QMFW.ml)) ./ 

(2~s SW(1 : 2" (p-pmax) : 2~pmax , 1 : 2~ (p-pmax) : 2~pmax) ; 

# 

# Computation of the point value filter related to TAUTILDE 

# 

TAUTW = ( (2~(ps) * PHIW* * P'HlW) .* SS)./ SIGMA; 

TAUTW 5B Per iodize (TAUTW, p) ; 

The subroutine Periodize is not described here, but it is a straight forward 
transcription of 43. 

III.4.2 Main Program 

Main Program Ellip 

[UX] = Ellip (FX , QMFW , FIW , QMFBIW , TAUTW ) 

#INPUT : 

#FX -> 2D array containing the sampling of the function 

# f; Size(FX) -> (2~p X 2~p); FX(i,j) - f (i/2~p, j/2~p) , 

# (i, j ) belong to fO, . . . ,2~p-l}~2. 

#FIW -> ID tab of data containing the interpolation filter 

# related to PHI_pO, size (FIW ) ->2~p ; 

#QMFW -> see PRECAL 

#QMFBIW ~> see PRECAL 
#TAUTW -> see PRECAL 

# 


quadrature mirror filters associated 

to the scaling functions; 

size (QMFBIW . PTILDEO)-> (2~p X 2~p) : 

QMFBIW. PTILDEl-> 2D tab containing the biorthogonal 
filters associated to the first wav- 
elet ; size (QMFBIW. PTILDE l)-> (2~p X 2 ~p). 

QMFBIW. PTILDE2-> same as QMFW. PTILDE 1 for the second 
wavelet (not computed QMFBIW . PTILDE2 = 
QMFW. PTILDE 1 transposed). 

QMFBIW. PTILDE3-> same as QMFW.PTILDE1 for the third 
wavelet . 
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# 

#OUTPUT : 

#UX ~> 2D array containing the sampling of the approxi- 

# -mation u_p; size(UX)->(2~p X 2~p), UX(i,j) = 

# u(i/2~p, j/2~p) , (i,j) belong to {0, . . . ,2~p-l}~2. 

# 

#TAUTW -> 2D array containing the values of the Fourier 

# transform of the scaling function TAUTILDE at level p 

# 

#TEMP0RARY DATA: 

#FW -> 2D tab containing the fft of FX; Size(FW) ■-> 

# (2~p X 2"p) ; 

#CPW -> 2D tab containing the fft of scaling coefficient of 

# FX; Size(CPW) -> (2~p X 2~p) ; 

#DJW -> Structure of 2D array containing the Fourier transform 

# of the wavelet coefficients; size(DJW) -> (2~p X 2~p); 
#CTILDEPW 

# -> same as CPU for UX; 

#UW -> 2D tab containing the fft of UX; 

# 

# 

# step 0 

# 

[FW] ■= Fast Fourier Transform(FX) 

# 

# step 1 

# 

[CPU] = FW.+FIW (Term by term multiplication) 

# 

# step 2 

# 

[DJW] = 2D Tensorial Tree Algorithm^ (CPW,QMFW) 

# 

# step 3 

# 

[DJW] = DJW.*(2~js) (Term by term multiplication) 

# 

# step 4 

# 

[CTILDEPW] = 2D Non Tensorial Tree Algorithm^ (D JW , QMFBI ) 

# 

# step 5 

# 
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[UW] = CTILDEPW . *TAUTW (Term by term multiplication) 

# 

[UX] = Inverse Fast Fourier Transform(UW) 

> 

III. 5 Storage and Complexity Analysis 

As the computation is clearly separated into precalculations and actual imple- 
mentation of the algorithm, we will also present the storage and complexity 
analyses separating the two parts. One should remember that the precalcu- 
lation is done once and for all while, as it will be the case in section IV, the 
algorithm can be applied iteratively. 

We will not discuss the complexity related to one-dimensional computations 
as well as the storage connected to one-dimensional arrays since both can be ne- 
glected in our bidimensional implementation. All the evaluations are performed 
for N = dim V p = 2 P x V, 

Storage 

Permanant storage (precalculations): The structures QMFBIW and TAUTW rep- 
resent four bidimensional arrays of size N. 

Temporary storage (actual algorithm): The storage related to bidimensional 
arrays can be reduced to one arrays of size N. 

Finally, the total memory used corresponds to five arrays of size N . 

Complexity analysis 

Precalculus: The computation of the four arrays in the structure QMFBIW is 
done in C x N. 

The computation of P V~ v is performed in C x N operations. The value of 

T p0 

C depends on the precision of the calculation. 

Main program: 

Fast Fourier Transform and Inverse Fast Fourier Transform in- 
volve C x N log(jV) operations. 

The complexity of the Term by term multiplication is N. 

Tree algorithm-D and Tree algorithm-1 are based on convolution and 
decimation operators. These procedures involve therefore C x N opera- 
tions. 

Therefore the total complexity is 0(Nlog(N)). 

In the following section we use these programs iteratively, to solve the 2D 
Burgers equation after reducing it to a cascade of elliptic problems. We would 


23 



like to emphasize moreover, that our approach can be also used to solve equa- 
tions involving homogeneous pseudo-differential operators. A characteristic ex- 
ample is \/-A u == / with periodi c bound ary conditions on [0, l] 2 . We have 
L = V— A and therefore <r(£) = + if* The most natural choice for S is 

S(£) = 2y sin 2 (£i/2) -fsin 2 (^2/2) and one easily checks that the hypotheses of 
theorem II.2 are then satisfied. The algorithms presented previously can be 
used (see Pj. Ponenti [16]). 

IV Numerical Application: Resolution of the 
2D Burgers equations 

In this section, we will use the periodized Battle- Lemarie’s multiresolution anal- 
ysis of splines of order m (see P.G. Lemarie [9]). The existence of collocation 
projectors related to the spline breakpoints requires splines of even order and 
the value m == 8 will be used in the applications. 

As described in J. Liandrat et al. [11], any parabolic equation of the type 

+ Lqu + G{u ) = 0 
u(Q,t) = u(l,t) 

< (44) 

u = uq for t = 0 

0 < t < Tmax , 0 < x < 1 

where L 0 is a differential operator of even order with positive symbol <7o(w), and 
G is generally a nonlinear function of u and its derivatives, can be numerically- 
approached using a classical finite difference time discretization scheme followed 
by a variational approximation of the resulting elliptic problems. We show now 
that the approach developed in the previous sections can be used efficiently to 
provide this approximation. 

Following J. Liandrat et al. [11] we first introduce a segmentation { t n }??~i 
of [O.Tmax] (i.e. a sequence {t n }n-i su °h that 0 = to < ti < ... < Im = Tmax) 
and now look for a sequence of functions of the x variable such 

that u^ n \x) is an approximation of u(x , t n ). 

With A t n = t n+1 — t n , 0 < n < M, and considering first (44) as an ordinary 
differential equation in time, a standard finite-difference discretization leads to 
the following iterative equation: 


£nU (n+1) = A t n , ..... A (45) 

where C n is a step forward operator that together with F is determined by 
the choice of the finite-difference approximation of the time-dependent ordinary 
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differential equation. We always assume that this approximation is at least semi 
implicit for the linear part Lq and explicit for the nonlinear part. Therefore 
3a > 0 such that C n = (I + aAt n Lo) where I stands for the identity operator. 
By hypothesis, cr 0 (w) > 0, Vo? and then C n has always a symbol bounded from 
below by 1. 

Hence, assuming that {«(""*),/ = 0,...,z} and {G(u( n ~ l \ l == are 

known, the resolution of (45) falls under the scope of paragraph III. 3 and the 
resolution of (44) can be therefore performed iteratively. 

The bidimensional Burgers equation writes, with u = (ui, u 2 ): 

-f Vu.u = uAu 
u(0,t) = u(M) 

< u == uq for t = 0 

0 < t < Tmax , 0 < x < 1. 

Choosing a constant step segmentation of [0, Tmax] (i.e., a segmentation 
such that EiAtf such that VO < n < M,t n = nAf), an implicit Crank-Nicholson 
time scheme (second order) for the linear term (z/Am), and an explicit second 
order Adams-Bashforth scheme for the nonlinear term (Vw.u) we get 

«( n+1 ) == + — 6 *(-— 

y 2 y 2 2 

and the solution can be written as 

u ( n+1 ) = fiC”+ 1 ) - -u(") , 


4 a 


with 

(47) 

To fall completely under the scope of paragraph III.3, one should be able 
to evaluate the point values of the nonlinear term of (47). We used the sim- 
plest method available that consists, as classically done in spectral methods 
(C. Canute et al. [3]), to “apply the nonlinear operator on the grid points”. 
More precisely, the approximation of we used is P7(Vw®.u®) = 

Cyv(C v v x (yu l ).Cyv{u 1 )) where x is a term by term multiplication of finite 
sequences and Cyv is the collocation projection introduced in section III.2. 

Then, for each time step nAt, the problem clearly belongs to the class of 
elliptic problems studied in section (III) with L = I — v^A and / == 2u.( n ) — 
£f(§Vtt( n ).u( n ) — The iterative form of the equation induces 

some modifications of the general scheme presented in III. 3 and we therefore 
provide the full scheme for an iteration of the Burgers approximation scheme: 


w( n+1 ) - (j - v | a) ^2u< rt * - i 
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0 . The inputs are the values of («( n )(« 7 p )), (^7 ~(J P )) and r(^p ))• 

1 . F^ n \J p ) = 2u^ n ' ) (J p ) — 6 t PV{Vu^ n \u^)(J p ) is computed as described 
above. 

2 . F p (x), the function of V p interpolating the values F n (J p ) is constructed 
as Fp = Hk^I p C P^pk > 

3 . F™ is then decomposed into the wavelet subspaces WJ , 0 < j < p — 1 and 

vy as 

F;= X) (I?,*a>®a+coo 

A6AJ- 1 

where coo = II ^(F"). 

4. 4" +1) then becomes 

4” +1) = E 2-> s {f;,^>^ + coo 

AeAj- 1 


where 0% = {V s L~ l i>x) V ■ 

5. 4 n+1) is then expanded in terms of the set of functions {j pk , fe € Ip} as 

4 n+1) = E • 

fee/. 

We also get ^u (n+1) = £ fc€/p c pk^ k 

= Hjfeg/y Cpfc^Ppfc' 

6. From the point values of t£(J p ), ^■r J f(J p ), and -§^r p (J p ) we compute 
(u( n+1 )(J p )) and its gradient using the corresponding point value filters. 

7. Finally using the values of ( u„(J p )) and of its gradient, we get the values 
of (u( n+1 )(Jp)) = (5(" +1 >(J p ) - uW(J p )) and its gradient, 

IV. 1 Storage and Complexity Analysis 

As described above, the numerical code implements the elliptic solver in an it- 
erative process. Since the time step At is constant, the characteristics of the 
elliptic solver does not depend on the time index n. Then, the solver precalcu- 
lations related to the whole parabolic problem are the same as the ones related 
to the elementary elliptic solver (see section III.5). This applies to the storage 
and to the complexity as well. 
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As it has been shown in the previous section however, extra work, not con- 
nected to the elliptic solver itself but to the computation of the right hand side 
term of the iterative equation (47) is required. This extra work is related to the 
storage of the fields at the different time steps involved in the three level time 
step Adams-Bashforth Crank- N icholson scheme and to the point value evalu- 
ation of the derivatives involved in the nonlinear part. Again, it can be split 
into permanent and temporary storage as well as in precalculation and main 
program extra work. 

Storage (in addition to the elliptic solver storage) 

Permanant storage (precalculations): One extra structure containing the 
point values ^r^(J p ) must be stored in one bidimensional array of size N; the 
structure containing the point values ^r p (J p ) is given by transposition of the 
previous one. 

Temporary storage ( actual algorithm): The two fields u^ n \ and 

can be handled using three arrays of size N. 

Complexity analysis (in addition to the elliptic solver complexity) 

Precalculation: The computation of PV d ~ v is performed in C x N opera- 

ax ~p * 

tions where, as in section III.5 the value of C depends on the precision of the 
calculation. 

Main program: The addition of complexity is related to steps (1), (5), (6) 
and (7). Since these steps involve convolutions and term by term products, the 
added complexity is again CN log N 

Finally, the total memory used is 7 arrays of size N- . 

The total complexity is Q(N log(N)) . 

Obviously, the total complexity of the whole resolution is M times the com- 
plexity of one time step resolution. 

IV. 2 Numerical Results 

Test case on an x 2 translation invariant problem 

The validation of our code has been performed on an x 2 translation invari- 
ant problem constructed using for the initial condition ( uo 11 uq 2 )( xi , x 2 ) == 
(sin(27r#i), 0). Indeed, with such an initial condition, the solution remains x 2 
translation invariant. 

For an easy comparison to the well documented paper of C. Basdevant et al, 
[1] we used v = 10" 2 /tt. 

As explained in C. Basdevant et al. [1], the pertinent quantities are 

Ou Ou 

ms = suptll— (z^^lloo = sup te [o,Tm a x) \ (0.5,*)| 
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and t ms defined such that 


|^(0.5,W)| = ms. 

Table 48 exhibits the numerical results obtain using various values for the 
time step At. The maximum time step numerically acceptable was At = 0.0075. 
In each case, the values of ms are computed by interpolation and the correspond- 
ing values of t ms are deduced. The comparison with the expected theoretical 
values (first column) shows that our method competes favorably with the ma- 
jority of the schemes presented in C. Basdevant et al. [1]. A complete study of 
the time step size dependence of the results connected to the stability analysis 
of the parabolic algorithm will be published later. 


Exact 


0.0005 

0.001 

0.0025 

0.005 

-304.0103 

ms 

-304.6308 

-305.727 

-309.4354 

-316.5454 

0.255237 

tms 

0.253 

0.252 

0.25 

0.245 


Test case on a first diagonal translation invariant problem 
Our second test case is performed on a first diagonal translation invariant prob- 
lem constructed using (u 0l , tio 3 ).(*ij*2) = (sin(27r(#i + ar 2 ),sin(27r(ici + # 2 )). 
Again, the solution can be compared to the reference solution of C, Basdevant 
et al. [1] thanks to a 45° rotation and to a time dilation of factor 2. However, 
according to our reference axes, it is obviously a fully bidimensional solution. 

Figures 1,2 and 3 show the isoline values of the numerical approximations 
computed with St — .001 at t = 0, t = 0,15 and t = 0.50. The first diagonal 
translation invariance is kept and we obtain the values ms — 249.0528 and t ms = 
0.123. The expected theoretical values are —304.0103 for ms and 0.1276185 for 
t m$ . This is not as good as before but one should note that the resolution in 
the direction perpendicular to the front axis is now half the one in our previous 
calculations. 

Since the ultimate application of all this work is the development of adaptive 
algorithms (i.e. the development of algorithms handling approximation spaces 
of reduced dimension adapted to the solution regularity (see for instance Pj. 
Ponenti [16]), we have estimated, at various times, the wavelet basis adapted 
to the approximation and defined as the lowest cardinal m = 8 spline wavelet 
basis preserving the L 2 norm of the approximation with a precision of 10~“ 6 . 
The columns of table (49) show for each scale 0 < j < 7 the number of wavelets 
selected in the adapted basis related to the approximated solution at various 
times. It appears that, compared to the full basis of Vs (last column), these 
bases have a drastically reduced cardinal (we defined the rate of compression 
as cardma ' c ^ n afj^ ^ ) even ^ the gradients of the solution fill up a large 
domain made of two complete lines of the plane (see figures 1, 2 and 3) 
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Figure 2: Approximated solution, t* — 0.15. 
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Figure 3: Approximated solution, ti = 0.50. 
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Conclusion 

In this paper, we have proposed an inversion scheme for elliptic problems 
based on biorthogonal wavelets. The approximation of elliptic problem solutions 
is constructed and leads to stable and fast numerical algorithms. 

Numerical tests related to the approximation of the parabolic Burgers equa- 
tions transformed into a cascade of elliptic problems are provided. 

The approximation scheme is based on convolution operators and can there- 
fore be theoretically used in the framework of adapted spaces of approximation. 
As mentioned however, the nice tensorial product structure that enforces nu- 
merical efficiency is fragile and is generally destroyed when applying the scheme 
directly. Other approximations for the step forward operator, should allow 
one to use efficiently this approximation in a general context of adapted multi- 
dimension spaces of approximation. 
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